####################################################
### Estimate Supply Model for Covered California ###
####################################################
	
##### Simulation parameters
	
	year = 2018;
	ns = 50;
	add_insurer_market_fe = true;
	add_av_interactions = true;
	inattention = true;
	network_flag = false;
	eq_robust_flag = false;
	control_function = true;
	random_parameters = true;
	keep_two_choices = false;
	separate_metal_flag = false;
	add_complex_interactions = true;
	nested = false;
	
	choice_seq_flag = false;	
	av_rs_flag = false;
	int_rs_variables = true;
	all_rs_variables = false;
	old_flag = false;
	hisp_flag = false;
	agegender_flag = false;
	if old_flag
		choice_seq_flag = false;
		separate_metal_flag = false;
	end
	
##### Load Data
	
	cd("C:/Users/eas24f/OneDrive - Florida State University/Research/Inertia RESTAT/Code")
	include("load.supply.estimation.data.jl")
	
##### Define covariates for regressions	
	
	all_covariates = convert(Array,param_estimates[!,:Variable]);
	if nested
		pop!(all_covariates);
	end
	product_covariates = intersect(all_covariates,
		["Premium","AV","HMO","Anthem","Blue_Shield","Kaiser","Health_Net"]);
	rating_area_variables = ["share_ra2","share_ra3","share_ra4","share_ra5","share_ra6","share_ra7",  
				"share_ra8","share_ra9","share_ra10","share_ra11","share_ra12","share_ra13", 
				"share_ra14","share_ra15","share_ra16","share_ra17","share_ra18","share_ra19"];
	claims_moments_covariates = cat(dims=1,["intercept","HMO","log_risk_score","trend",
		"Anthem","Blue_Shield","Kaiser","Health_Net"],rating_area_variables);
	
	if all_rs_variables
		risk_score_share_variables = cat(dims=1,["share_0to34","share_35to54","share_250to400","share_gt400",
			"share_male","share_family","share_Asian","share_Black","share_Hispanic","share_Other_Race"]);
	elseif int_rs_variables
		risk_score_share_variables = cat(dims=1,["share_0to34","share_male","share_family","share_minority"]);
	elseif hisp_flag
		risk_score_share_variables = cat(dims=1,["share_0to34","share_male","share_family","share_Hispanic"]);
	elseif agegender_flag
		risk_score_share_variables = cat(dims=1,["share_0to34","share_male"]);
	elseif old_flag
		risk_score_share_variables = cat(dims=1,["share_18to54","share_Hispanic"]);
	end
	
	if av_rs_flag
		risk_score_covariates = cat(dims=1,"intercept","AV",risk_score_share_variables);
	else
		risk_score_covariates = cat(dims=1,"intercept","Silver","Gold","Platinum",risk_score_share_variables);
	end
	
	if inattention
		premium_shock_index = findall(all_covariates .== "premium_shock")[1];
		search_covariates = all_covariates[premium_shock_index:length(all_covariates)];
		all_covariates = all_covariates[1:(premium_shock_index-1)];
		
		beta_search = betap[premium_shock_index:(length(betap))];
		min_beta = betap[1:(premium_shock_index - 1)];
		
		control_function = false;
	end
	
	all_mkts = collect(2:1:19);
	premium_indices = findall(convert(Array,param_estimates[!,:Premium_index]) .== 1);	
		
		
##### Extract sample - only what the firm knows in year "year"

	keep_households = findall(convert(Array,households[!,:year]) .<= year);
	if eq_robust_flag
		keep_households = intersect(keep_households,findall(convert(Array,households[!,:year]) .>= 2016));
	end
	
	if keep_two_choices
		years_in_panel = combine(groupby(households[keep_households,:],:household_id), :household_id => length)[!,:household_id_length];
		time_in_panel = zeros(size(households)[1]);
		unique_hseq_ids = unique(households[keep_households,:household_id]);
		for s in 1:length(unique_hseq_ids)
			time_in_panel[findall(households[!,:household_id] .== unique_hseq_ids[s])] .= years_in_panel[s];
		end
		
		# Delete households only in panel for 1 year
		keep_households = intersect(keep_households,findall(time_in_panel .> 1));
			
		# Delete households' choices in years 3 and above
		first_year_in_panel = combine(groupby(households[keep_households,:],:household_id), :year => minimum)[!,:year_minimum];
		in_first_two_years = zeros(size(households)[1]);
		unique_hseq_ids = unique(households[keep_households,:household_id]);
		for s in 1:length(unique_hseq_ids)
			in_first_two_years[intersect(findall(households[!,:household_id] .== unique_hseq_ids[s]),findall(households[!,:year] .<= first_year_in_panel[s] + 1))] .= 1;
		end
		keep_households = intersect(keep_households,findall(in_first_two_years .== 1));
		
		# Some Households are in sample for two non-consecutive years - delete them
		years_in_panel = combine(groupby(households[keep_households,:],:household_id), :household_id => length)[!,:household_id_length];
		time_in_panel = zeros(size(households)[1]);
		unique_hseq_ids = unique(households[keep_households,:household_id]);
		for s in 1:length(unique_hseq_ids)
			time_in_panel[findall(households[!,:household_id] .== unique_hseq_ids[s])] .= years_in_panel[s];
		end
		keep_households = intersect(keep_households,findall(time_in_panel .> 1));
	end	
	
	new_household_numbers = zeros(Int64,size(households)[1]);
	new_household_numbers[keep_households] = collect(1:1:length(keep_households));
	data[!,:household_number] = new_household_numbers[convert(Array,data[!,:household_number])];
	
	data = data[data[!,:household_number] .> 0,:];
	households = households[keep_households,:];
	if eq_robust_flag
		data = data[findall(convert(Array,data[!,:year]) .>= 2016),:];
		households = households[findall(convert(Array,households[!,:year]) .>= 2016),:];
	end
	
	uninsured_plans = convert(Array,data[!,:uninsured_plan]);
	household_numbers_data = convert(Array,data[!,:household_number]);	
	household_numbers = household_numbers_data[uninsured_plans .== 0];	
	choices = convert(Array,data[!,:choice]);
	
	household_weights = data[uninsured_plans .== 0,:weight];
	
	##### Remove any plans not available in sample
	
	data_pmt_names = string.(data[uninsured_plans .== 0,:plan_name],
								data[uninsured_plans .== 0,:year],
								households[household_numbers,:rating_area]);
	data_pt_names = string.(data[uninsured_plans .== 0,:plan_name],
								data[uninsured_plans .== 0,:year]);
	
	plan_names = sort(unique(data[uninsured_plans .== 0,:plan_name]));
	plan_year_names = sort(unique(data_pt_names));
	plan_pmt_names = sort(unique(data_pmt_names));
	
	plan_indicator = trues(size(plans)[1]);
	for i in eachindex(plan_indicator)
		plan_indicator[i] = in(plans[i,:pt_name],plan_year_names);
	end
	plans = plans[findall(plan_indicator),:];
	sort!(plans);
	
	plan_indicator = trues(size(plans_pmt)[1]);
	plans_pmt[!,:pmt_name] = string.(plans_pmt[!,:plan_name],plans_pmt[!,:year],plans_pmt[!,:rating_area]);
	for i in eachindex(plan_indicator)
		plan_indicator[i] = in(plans_pmt[i,:pmt_name],plan_pmt_names);
	end
	plans_pmt = plans_pmt[findall(plan_indicator),:];
	sort!(plans_pmt);
	
	insurer_vector = convert(Array,plans[!,:Insurer]);
	insurers = unique(insurer_vector);
	
	household_pmt_indices = zeros(Int64,length(data_pmt_names));
	own_plan_year_matrix = zeros(Int64,length(plan_pmt_names),length(plan_pmt_names));
	plan_year_matrix = zeros(Int64,length(plan_pmt_names),length(plan_pmt_names));
	firm_plan_matrix = zeros(Int64,length(insurers),length(plan_pmt_names));
	for j in 1:length(plan_pmt_names)	
		firm = plans_pmt[j,:Insurer];
		yearp = plans_pmt[j,:year];
		mkt = plans_pmt[j,:rating_area];
		own_plan_year_indices = intersect(findall(plans_pmt[!,:Insurer] .== firm),
			findall(plans_pmt[!,:year] .== yearp));
		plan_year_indices = findall(plans_pmt[!,:year] .== yearp);
		plan_indices = findall(data_pmt_names .== plan_pmt_names[j]);
			
		firm_plan_matrix[findall(insurers .== firm)[1],j] = 1;	
		own_plan_year_matrix[own_plan_year_indices,j] .= 1;
		plan_year_matrix[plan_year_indices,j] .= 1; 
		household_pmt_indices[plan_indices] .= j;
	end
	
	household_pt_indices = zeros(Int64,length(data_pt_names));
	for j in 1:length(plan_year_names)
		plan_indices = findall(data_pt_names .== plan_year_names[j]);
		household_pt_indices[plan_indices] .= j;
	end
	
	year_plan_names = string.(plans_pmt[!,:plan_name],plans_pmt[!,:year]);
	year_plan_indices = zeros(Int64,size(plans_pmt)[1]);
	for j in 1:size(plans)[1]
		year_plan_indices[intersect(findall(plans_pmt[!,:plan_name] .== plans[j,:plan_name]),
			findall(plans_pmt[!,:year] .== plans[j,:year]))] .= j;
	end
	
	year_plan_moment_names = string.(moments_data[!,:plan_name],moments_data[!,:year]);
	year_plan_moment_indices = zeros(Int64,size(moments_data)[1]);
	for j in 1:size(plans)[1]
		year_plan_moment_indices[intersect(findall(moments_data[!,:plan_name] .== plans[j,:plan_name]),
			findall(moments_data[!,:year] .== plans[j,:year]))] .= j;
	end
	
	plan_foc_moment_cond_names = string.(plans_pmt[!,:plan_name],plans_pmt[!,:year]);
	foc_moment_cond_names = sort(unique(plan_foc_moment_cond_names));
	foc_moment_indicator = zeros(Int64,length(plan_pmt_names),length(foc_moment_cond_names));
	for m in 1:length(foc_moment_cond_names)
		plan_indices = findall(plan_foc_moment_cond_names .== foc_moment_cond_names[m]);
		foc_moment_indicator[plan_indices,m] .= 1;
	end
	
	if random_parameters 
		
		if inattention
			random_variables = ["Anthem","Blue_Shield","Kaiser","Health_Net"];
			random_parameter_product_indices = findall(in(random_variables),product_covariates);
			random_product_covariates = string.(random_variables,"_random");
		elseif add_av_interactions
			random_variables = ["intercept","Anthem","Blue_Shield","Kaiser","Health_Net"];
			random_parameter_product_indices = cat(dims=1,findall(all_covariates .== "intercept")[1],findall(in(random_variables),product_covariates));
			random_product_covariates = string.(random_variables,"_random");
		else
			random_variables = ["intercept","Anthem","Blue_Shield","Kaiser","Health_Net","AV"];
			random_parameter_product_indices = cat(dims=1,findall(all_covariates .== "intercept")[1],findall(in(random_variables),product_covariates));
			random_product_covariates = string.(random_variables,"_random");
		end
		
		
		random_parameter_indices = findall(in(random_product_covariates),all_covariates);
		if length(random_parameter_indices) == 1
			random_parameter_indices = random_parameter_indices[1];
		end
		
		# Assign a household reference number
		household_ref_numbers = zeros(maximum(households[!,:household_id]));
		household_ref_numbers[unique(households[!,:household_id])] = collect(1:1:length(unique(households[!,:household_id])));
		households[!,:ref_number] = Int64.(household_ref_numbers[households[!,:household_id]]);
		
		V = zeros(size(data)[1],length(random_variables),ns);
		for k in 1:length(random_variables)
			Random.seed!(random_parameter_product_indices[k]);
			Vi = randn(length(unique(households[!,:ref_number])),ns);
			Vi = Vi[households[!,:ref_number],:];
			V[:,k,:] = Vi[household_numbers_data,:];
		end
	end
	
	moments_data[!,:Insurer] .= "Small_Insurer";
	moments_data[moments_data[!,:Anthem] .== 1,:Insurer] .= "Anthem";
	moments_data[moments_data[!,:Blue_Shield] .== 1,:Insurer] .= "Blue_Shield";
	moments_data[moments_data[!,:Health_Net] .== 1,:Insurer] .= "Health_Net";
	moments_data[moments_data[!,:Kaiser] .== 1,:Insurer] .= "Kaiser";
	
	plan_ra_factors = ones(length(plan_pmt_names)) * 0.6;
	plan_ra_factors[plans_pmt[!,:AV] .== 0.7] .= 0.7 * 1.03;
	plan_ra_factors[plans_pmt[!,:AV] .== 0.9] .= 0.9 * 1.15;
	plan_ra_factors[plans_pmt[!,:AV] .== 0.8] .= 0.8 * 1.08;
	
	nu_for_riskadj = 0.86;
	
	metal_matrix = cat(dims=2,(plans_pmt[!,:AV] .== 0.6) .* 1,(plans_pmt[!,:AV] .== 0.7) .* 1,(plans_pmt[!,:AV] .== 0.8) .* 1,(plans_pmt[!,:AV] .== 0.9) .* 1);
	
##### Create Data Objects
	# hdata object: contains household data (things that don't change)
	# hvar object: contains household variables (things that change when we run counterfactuals)
	# pmt object: contains plan-market-year variables
	
	# X object: demand covariates
	# X_claims_data: claims covariates (data to use for estimation - not all plans are present b/c of missing data)
	# X_CLMS: claims covariates (data containing all predictions for all plans in years 2014 to t-1)
	# X_RS_data: risk score covariates (data to use for estimation - not all plans are present b/c of missing data)
	# X_RS: risk score covariates (data containing all predictions for all plans)
	# X_T_data: trend covariates (data to use for estimation - not all plans are present b/c of missing data)
	# X_T: trend covariates (data containing all predictions for all plans in year t)
	
	# Field names
		
		hdata_fields = ["penalty","weight","pricing_factor","rating_factor","age_factor","ra_factor","GCF","alpha",
						"real_prob","real_ex_prob","real_instrumented_prob","risk_score","members","choice","real_slc_plan","real_unsubs_premium"];
		hvar_fields = ["subsidized_premium","subsidy","slc_plan","prob","ex_prob","risk_score"];
		pmt_fields = ["risk_score","ra_factor","real_demand","plan_pricing_weight","average_age_factor",
						"total_premiums","reins_factor","variable_admin","fixed_admin","base_premium"];
		partial_names = ["partial_plan_demand_partial_price","partial_firmrevenue_partial_price",
			"partial_firmtransfer_partial_price","partial_firmadmin_partial_price"];
		
		X = form_X_matrix(data,households,plans,product_covariates); 
		if inattention
			Xsearch = form_Xsearch_matrix(search_covariates,households);
		end		
		X_claims_data = form_X_claims_matrix(moments_data,claims_moments_covariates);
		
		y_claims_data = convert(Array,moments_data[moments_data[!,:year] .<= year,:log_cost]);
		y_RS_data = convert(Array,moments_data[moments_data[!,:year] .<= year,:log_risk_score]);
		
		weights_data = convert(Array,moments_data[moments_data[!,:year] .<= year,:EXP_MM]);
		hh_share_matrix = update_dem_share_matrix(households);
		ac_weights = rescale_weights(moments_data,weights_data,2014,year);
				
		hdata = initialize_household_data(X,hdata_fields,data,plans,household_numbers,min_beta);
		pmt = initialize_pmt_data(pmt_fields,plans_pmt,hdata); 
		X_RS = form_X_RS_matrix(plans_pmt,hdata,household_numbers,risk_score_covariates);
		X_CLMS = form_X_CLMS_matrix(plans_pmt,hdata,claims_moments_covariates,rs_estimates[!,2]);
		
		# Weights for moment contributions
		# Weights for plans that map to same proj_name
		moment_cont_weights = zeros(size(plans_pmt)[1]);
		plan_pmt_demand = pmt[:,findall(pmt_fields .== "real_demand")[1]];
		for j in 1:size(moments_data)[1]
			pmt_indices = findall(plans_pmt[!,:rate_name] .== moments_data[j,:pmt_name]);
			moment_cont_weights[pmt_indices] = ac_weights[j] .* 
				(plan_pmt_demand[pmt_indices] ./ sum(plan_pmt_demand[pmt_indices]));
		end
		
#### Now we do GMM to estimate the cost parameters
	# Could estimate both demand and cost jointly here, but this would be very computationally intensive (take demand parameters from MLE)
	# Moments:
		# Cost moments
		# Risk score moments
		# Predicted cost moments
		# FOCs
	
	# First, let's take the risk score parameters as given and estimate the cost parameters
	
		num_moment_conditions = length(claims_moments_covariates) + length(foc_moment_cond_names);
		
		#theta_CLMS_initial = (CSV.read(string("claims_initial_",year-1,".csv"),header=true)[!,2])[1:length(claims_moments_covariates)];
		initial_claims_estimates = CSV.read(string("claims_initial_",2017,".csv"),DataFrame,header=true);
		theta_CLMS_initial = convert(Array,initial_claims_estimates[!,2]);
		if separate_metal_flag
			theta_RS_initial_bronze = rs_estimates_bronze[!,2];
			theta_RS_initial_silver = rs_estimates_silver[!,2];
			theta_RS_initial_gold = rs_estimates_gold[!,2];
			theta_RS_initial_platinum = rs_estimates_platinum[!,2];
			theta_RS_initial = cat(dims=2,theta_RS_initial_bronze,theta_RS_initial_silver,theta_RS_initial_gold,theta_RS_initial_platinum);
		else
			theta_RS_initial = rs_estimates[!,2];
		end
		
		theta_initial = theta_CLMS_initial;
		theta_RS_share_indices = findall(in(risk_score_share_variables),risk_score_covariates);
		
		println("Calculate partial derivatives")
		if nested
			partials = calculate_key_partials(theta_RS_initial,pmt,hdata,X_RS,min_beta[length(min_beta)]);
		else
			partials = calculate_key_partials(theta_RS_initial,pmt,hdata,X_RS,1);
		end
		foc_reins_factor = pmt[:,findall(pmt_fields .== "reins_factor")[1]];
		
		# Extract from partials object
		jacobian = partials[1:length(plan_pmt_names),1:length(plan_pmt_names)];
		rs_partial_matrix = partials[1:length(plan_pmt_names),(length(plan_pmt_names) + 1):(2*length(plan_pmt_names))];
		partial_plan_demand_partial_price = partials[:,2*length(plan_pmt_names) + findall(partial_names .== "partial_plan_demand_partial_price")[1]];
		partial_firmdemand_partial_price = vec(sum(dims=1,jacobian .* own_plan_year_matrix));
		partial_firmrevenue_partial_price = partials[:,2*length(plan_pmt_names) + findall(partial_names .== "partial_firmrevenue_partial_price")[1]];
		partial_firmtransfer_partial_price = partials[:,2*length(plan_pmt_names) + findall(partial_names .== "partial_firmtransfer_partial_price")[1]];
		partial_firmadmin_partial_price = partials[:,2*length(plan_pmt_names) + findall(partial_names .== "partial_firmadmin_partial_price")[1]];	
		
		partial_firmclaims_partial_price_focs = (1 ./ partial_plan_demand_partial_price) .*
			(partial_firmrevenue_partial_price .+ partial_firmtransfer_partial_price .- partial_firmadmin_partial_price) ./ (1 .- foc_reins_factor);
	
		plan_risk_scores = exp.(X_RS * theta_RS_initial);
		plan_demands = pmt[:,findall(pmt_fields .== "real_demand")[1]];
		
		
		# Step 1 - Initial guess for theta_CLMS and theta_RS (I just did OLS to get this)
		W = Matrix{Float64}(I,num_moment_conditions,num_moment_conditions);
		println("Start Step 1 Estimation")
		step1_results = optimize(compute_gmm_objective,theta_initial,iterations=1000000,show_trace=true);
		theta_step1 = Optim.minimizer(step1_results);
		if inattention
			CSV.write(string("step1_results_inattention.csv"),
				DataFrame(cat(dims=2,cat(dims=1,claims_moments_covariates),theta_step1),:auto));
		elseif network_flag
			CSV.write(string("step1_results_net.csv"),
				DataFrame(cat(dims=2,cat(dims=1,claims_moments_covariates),theta_step1),:auto));
		elseif eq_robust_flag
			CSV.write(string("step1_results_eqrobust.csv"),
				DataFrame(cat(dims=2,cat(dims=1,claims_moments_covariates),theta_step1),:auto));
		elseif (old_flag & av_rs_flag)
			CSV.write(string("step1_results_old_av.csv"),
				DataFrame(cat(dims=2,cat(dims=1,claims_moments_covariates),theta_step1),:auto));
		elseif old_flag
			CSV.write(string("step1_results_old.csv"),
				DataFrame(cat(dims=2,cat(dims=1,claims_moments_covariates),theta_step1),:auto));
		elseif (int_rs_variables & av_rs_flag)
			CSV.write(string("step1_results_int_av.csv"),
				DataFrame(cat(dims=2,cat(dims=1,claims_moments_covariates),theta_step1),:auto));
		elseif (choice_seq_flag & int_rs_variables & add_av_interactions)
			CSV.write(string("step1_results_int_seq_avint.csv"),
				DataFrame(cat(dims=2,cat(dims=1,claims_moments_covariates),theta_step1),:auto));
		elseif (choice_seq_flag & int_rs_variables)
			CSV.write(string("step1_results_int_seq.csv"),
				DataFrame(cat(dims=2,cat(dims=1,claims_moments_covariates),theta_step1),:auto));
		elseif (int_rs_variables & add_av_interactions)
			CSV.write(string("step1_results_int_avint.csv"),
				DataFrame(cat(dims=2,cat(dims=1,claims_moments_covariates),theta_step1),:auto));
		elseif int_rs_variables
			CSV.write(string("step1_results_int.csv"),
				DataFrame(cat(dims=2,cat(dims=1,claims_moments_covariates),theta_step1),:auto));
		elseif hisp_flag
			CSV.write(string("step1_results_hisp.csv"),
				DataFrame(cat(dims=2,cat(dims=1,claims_moments_covariates),theta_step1),:auto));
		elseif agegender_flag
			CSV.write(string("step1_results_agegender.csv"),
				DataFrame(cat(dims=2,cat(dims=1,claims_moments_covariates),theta_step1),:auto));
		elseif (all_rs_variables & choice_seq_flag & add_av_interactions)
			CSV.write(string("step1_results_full_seq_avint.csv"),
				DataFrame(cat(dims=2,cat(dims=1,claims_moments_covariates),theta_step1),:auto));
		elseif (all_rs_variables & add_av_interactions)
			CSV.write(string("step1_results_full_avint.csv"),
				DataFrame(cat(dims=2,cat(dims=1,claims_moments_covariates),theta_step1),:auto));
		elseif (all_rs_variables & choice_seq_flag & av_rs_flag)
			CSV.write(string("step1_results_full_seq.csv"),
				DataFrame(cat(dims=2,cat(dims=1,claims_moments_covariates),theta_step1),:auto));
		elseif (all_rs_variables & av_rs_flag)
			CSV.write(string("step1_results_full_av.csv"),
				DataFrame(cat(dims=2,cat(dims=1,claims_moments_covariates),theta_step1),:auto));
		elseif all_rs_variables
			CSV.write(string("step1_results_full.csv"),
				DataFrame(cat(dims=2,cat(dims=1,claims_moments_covariates),theta_step1),:auto));
		end
		
		# Step 2
		println("Compute Optimal Weight Matrix")
		W = compute_optimal_weight_matrix(theta_step1);
		println("Start Step 2 Estimation")
		step2_results = optimize(compute_gmm_objective,theta_step1,iterations=60000,show_trace=true) ;
		theta_opt = Optim.minimizer(step2_results);
		se = compute_standard_errors(theta_opt);
		pval = 2 * (1 .- cdf.(Normal(),abs.(theta_opt ./ se)));
		if inattention
			CSV.write(string("claims_estimates_int_avint_inattention.csv"),
					DataFrame(cat(dims=2,cat(dims=1,claims_moments_covariates),theta_opt,se,pval),:auto));
		elseif network_flag
			CSV.write(string("claims_estimates_int_avint_net.csv"),
					DataFrame(cat(dims=2,cat(dims=1,claims_moments_covariates),theta_opt,se,pval),:auto));
		elseif eq_robust_flag
			CSV.write(string("claims_estimates_int_avint_eqrobust.csv"),
					DataFrame(cat(dims=2,cat(dims=1,claims_moments_covariates),theta_opt,se,pval),:auto));
		elseif (old_flag & av_rs_flag)
			CSV.write(string("claims_estimates_old_av.csv"),
					DataFrame(cat(dims=2,cat(dims=1,claims_moments_covariates),theta_opt,se,pval),:auto));
		elseif old_flag
			CSV.write(string("claims_estimates_old.csv"),
					DataFrame(cat(dims=2,cat(dims=1,claims_moments_covariates),theta_opt,se,pval),:auto));
		elseif (int_rs_variables & av_rs_flag)
			CSV.write(string("claims_estimates_int_av.csv"),
					DataFrame(cat(dims=2,cat(dims=1,claims_moments_covariates),theta_opt,se,pval),:auto));
		elseif (choice_seq_flag & int_rs_variables & add_av_interactions)
			CSV.write(string("claims_estimates_int_seq_avint.csv"),
					DataFrame(cat(dims=2,cat(dims=1,claims_moments_covariates),theta_opt,se,pval),:auto));
		elseif (choice_seq_flag & int_rs_variables)
			CSV.write(string("claims_estimates_int_seq.csv"),
					DataFrame(cat(dims=2,cat(dims=1,claims_moments_covariates),theta_opt,se,pval),:auto));
		elseif (int_rs_variables & add_av_interactions)
			CSV.write(string("claims_estimates_int_avint.csv"),
					DataFrame(cat(dims=2,cat(dims=1,claims_moments_covariates),theta_opt,se,pval),:auto));
		elseif int_rs_variables
			CSV.write(string("claims_estimates_int.csv"),
					DataFrame(cat(dims=2,cat(dims=1,claims_moments_covariates),theta_opt,se,pval),:auto));
		elseif hisp_flag
			CSV.write(string("claims_estimates_hisp.csv"),
					DataFrame(cat(dims=2,cat(dims=1,claims_moments_covariates),theta_opt,se,pval),:auto));
		elseif agegender_flag
			CSV.write(string("claims_estimates_agegender.csv"),
					DataFrame(cat(dims=2,cat(dims=1,claims_moments_covariates),theta_opt,se,pval),:auto));
		elseif (all_rs_variables & choice_seq_flag & add_av_interactions)
			CSV.write(string("claims_estimates_full_seq_avint.csv"),
					DataFrame(cat(dims=2,cat(dims=1,claims_moments_covariates),theta_opt,se,pval),:auto));
		elseif (all_rs_variables & add_av_interactions)
			CSV.write(string("claims_estimates_full_avint.csv"),
					DataFrame(cat(dims=2,cat(dims=1,claims_moments_covariates),theta_opt,se,pval),:auto));
		elseif (all_rs_variables & choice_seq_flag & av_rs_flag)
			CSV.write(string("claims_estimates_full_seq.csv"),
					DataFrame(cat(dims=2,cat(dims=1,claims_moments_covariates),theta_opt,se,pval),:auto));
		elseif (all_rs_variables & av_rs_flag)
			CSV.write(string("claims_estimates_full_av.csv"),
					DataFrame(cat(dims=2,cat(dims=1,claims_moments_covariates),theta_opt,se,pval),:auto));
		elseif all_rs_variables
			CSV.write(string("claims_estimates_full.csv"),
					DataFrame(cat(dims=2,cat(dims=1,claims_moments_covariates),theta_opt,se,pval),:auto));
		end	
		
		